


Institutional Archive of the Naval Postgraduate School 





Calhoun: The NPS Institutional Archive 
DSpace Repository 


Theses and Dissertations 1. Thesis and Dissertation Collection, all items 


1983-09 


SIMTBED a graphical test bed for analyzing 
and reporting the results of a statistical 
Simulation experiment 


Drueg, Hans-Walter 


Monterey, California. Naval Postgraduate School 
http://ndl.handle.net/10945/19984 
Copyright is reserved by the copyright owner 


Downloaded from NPS Archive: Calhoun 


Calhoun is the Naval Postgraduate School's public access digital repository for 
| (8 D U DLEY research materials and institutional publications created by the NPS community. 
«ist Ser Calhoun is named for Professor of Mathematics Guy K. Calhoun, NPS's first 


NY KNOX appointed — and published -- scholarly author. 

ies) LIBRARY Dudley Knox Library / Naval Postgraduate School 

411 Dyer Road / 1 University Circle 
Monterey, California USA 93943 





http://www.nps.edu/library 


] < rz 
. eet ete Le) 


tude Sonat Bove 


OF ea Le teh ed 
Ct et ne Lt et ee 
Ode pe ead Ded ed 


pom pyey yer fee yy mete 


ss ieaaeenl ane ioe oar onal 
ee tL le eS Pe ee) 
ie Lee ee eee a 


PO A eet ae eat et ae 
teed 
oii cry 


eo] te eed +g oe.) 
Sl are as DS Diivtheamtiet 14 
ae eee) ST) 


Cede ca te Aad ee Bed 


ial pape se hed gee ate Attach w= 


<r a ey ed eee 
een kee el es el ee eh i | 





CF Fehgh = eho oth"a. 
Ft 





Le 

; it Ma we 
ery Say a A tA Peper beret 
Lee Soe 
LEP LL OP 
7 pr o pi Spobs 

eee thd tie) 

4D aD 
Sat eens 


F) 





oe Se yy ee 
cr ene 





Pie 
Be ye 





| Put pilndatlonod-e Supls whedeted ott 


Tih a adichepatiel 1th ta 1 ee ae 
let eabepee prs 





LN ae stent ete ay A Kiet pe Le ee 
ore Dy 





tek anienaml Pun af tide cmon. 
ee Lee eee 
eh ae ae pe ee 
is ST. ee : 


5 
ee ae 
, tit Dandie. hatte mana 

a hee ee ee ee ee ee ee 


~~ 


ee 








Ce) 





cw ‘c- vo : 
Ce ee ee a er) ON Ee ere ve 





ae 
ys Oe ee ee 





EPEC ae Os 86 8 
CU PRs ies Ld 


ee See ee ee ot 5 
satin FRE So leh eRe oe ee ee ee oie 
PD cade ad a Le ee ee et ee eee OE atl = rts >». 
ee ie We le pi a eee a, 
ery Sates a. 
| ceseiahemics Skt Knee anne «tte eM a? sd 
ated Mig Yor aig aa sorter eelacties, teeta, 


ayy oP ve fe . 
Ls Oe Sat ee, Me en Se ee ee 
A 2 aah preperer terre ‘tee 
Leland Tah ee Vee ve 
2) Vy antes nq ee ee 
to Ce ee c 
Ce nt op hs eee oe, 


OPT rye sae 


eS ee ee 





| 


oe Le Tar semana spitesensl gs 
a ee ee Oe ee ree 
te a OR hedenane spud eo. 
ee ee a ey 


eres Se Rn Year ere = Aer: 4 
miami i ee a ee ee ee 
" bepderaned eG ecede th ok ate 


aa a ee 
pd edb note ve aerg porte 
eo vt ne ee 

map isiealeko: ponaa ae Se ae ee 


Sr rae ee 1 CTO Pe ee ee eee 
Siege ee es ana ee Nata a eRe Ye ee nt ee 7 “i . 


eee Sree f) CUAL Le caetleotes Me 
a - * 


ras 


= Ps OS PT 
Pa Oe artnet te Da hte ead een 
ee k 
oo. . 


ee bos 
a oe, ee 
Co pees are Ca ys 2 
pera re Peale 
Sed tS cag 
Saat) Te 
Seen Tenia eS 
P ai Ss aN 
ta J rn ee a ake 
Se a eas | 
Ce tins Sb, toes fheed ee ee 
ae Be, ae a Wem eee ve es ee 
Pal = - a Tg a be A ae ee 
\ a a ONG ate ee eee Te 
4 a we. ae eer 
cry 1° ee tee to ere ed 
C-— - wr] iF) o 
ca a ani aa es 
ca 7 wad x Le 
Filed be ee ee ah 
= sYhomethe af A Ce ee a) 
a 7 y a Chea Cal 
Peony as ee ley re 
Ate = ee 
a bY Cea te ne 
a. route in’ * Pe Oe eee ee et 
Fo a" i a 
a a F Ly 
I . rE 
Fumi Se vw ay 
Lick a rs) | ae Ce 
a ar ee ee eo) (aw. 
a fares . 
a ” a Ca 1 
= er a 
SA en SS L SEPA & 8 
: 4 J UL i C7 Se a 
noe si pees Se OAM Peis “Rae 
- Pt ad 
i « iy owe Cs 
Aaa 2 vy ov 
Ts . _ oy 
J ta we s reo bs 1 ee i 
~ ry 7 ae 
Ae ne es aes 
by 5 Ne 
ay tia ohne ot % 


is H 
ae) a 4 
Awsee & ite Fy Py 
<& hs, “ a ek 
ey ‘ 
eae b a aN 
. he oe Pt er ae 
aaah ts ae Ph 
> r 
ca Ld 2 a Sy 
Ce oY ae! 8 ee 
ry ye 
r ie a ok ra 
aye *alh, F CO ee 
LJ Ps mY od 
es 
> 
Cs A eke 
a 
oo] wy, we ee 
oo “ iret 
rn 
" Tr @« 
ron Cs fa al w 
Cd whe 
r - 
Ce) » 
ry bd ° ca 
ot 
5 ad e 
Sy 
ts - s 
Owe 
. 
7 & a 
os 
eo 
o C3 EN 
a me 
Sy 
eS 
rd 
vy a 
ed 


Ad 





Ce eo alee See SOLS to ; 
de ea Lien 





wre ro re c 
er yee evn : ah Meter Te Maes Peete ae ie ae 
ry ee ee k af oh arth tit Pte) 





Le et Pe eS 






Lo lee te Pek 
Ya aye Le ee ee 


See ere ay tee ere Shy 
ee ee ee te Y 

bo) ohenaetate ee lt 

A deen te et ee ee ee 








Ne ee ee yD | Re et ee te hn et 


Delt die anhenuttre Mota Meetalh tee beeen mt al 






ee er ee eer 


. Se aie ath teeth, Gn TE tte ed. 
ea Me rt thet de oe ee ee ae ee CS 
Se ee ee ee ee 
ON tlhe et te ee ee Oe ed 
CT De eh ee ee 


Oe ee he ee 
ein XW wen) des eT ee te tar ate op de ° 
. ry Ch a ee 


et ee 

















NAVAL POSTGRADUATE SCHOOL 


Monterey, California 





SIMTBED a Grapnical Test Bed for Analyzing 
and Revorting the Results of a 
Statisticai Simulation Experiment 







by 


Hans-Walter Crueg 






September 1983 





UNCLASSIFIED 


ee ap a 
SECURITY CLASSIFICATION OF THIS PAGE (When Data Entered) 


REPORT DOCUMENTATION PAGE SEP One ore rou 


- REPORT NUMBER 2. GOVT ACCESSION NO. 3. RECIPIENT'S CATALOG NUMBER 


4. TITLE (and Subtitie) 5S. TYPE OF REPORT & PERIOD COVERED 


SIMTBED a Graphical Test Bed for We eet Dh ass - 
Analyzing and Reporting the Results cane 1933 


Experiment 
. AUTHOR(a) 8. CONTRACT OR GRANT NUMBER(a) 
























Hans-Walter Drueg 









10. PROGRAM ELEMENT. PROJECT, TASK 
AREA & WORK UNIT NUMBERS 





- PERFORMING ORGANIZATION NAME ANO AOORESS 





Naval Postgraduate School 
Monterey, California 93943 









12. REPORT OATE 


September 1983 


» CONTROLLING OFFICE NAME ANO ADORESS 


Naval Postgraduate School 









13. NUMBER OF PAGES 













Monterey, California 93943 158 
. MONITORING AGENCY NAME & AOORESS(If dilferent from Controlling Office) 1S. SECURITY CLASS. (of thie report) 
Unclassified 










OECLASSIFICATION/ OOWNGRADOING 
SCHEOULE 


\Se. 
. OISTRIBUTION STATEMENT (of thie Report) 


Approved for public release; distribution unlimited. 






. BISTRIBUTION STATEMENT (of the ebetraci entered in Block 20, if different from Report) 









- SUPPLEMENTARY NOTES 





» KEY WOROS (Continue on reveree elde if neceeeary and identity by block number) 


SIMTBED Gamma astm HouLION 
Seatistical Simulation Standard Deviation of 
Graphics Corne latromecoer f1ClLene 
Regression-Adjusted Graphics 
Serial Correlation Coefficient 
20. ABSTRACT (Coniinue on reveree aide i{ neceseary and identify by dlock number) 


| A graphical test bed in which the results of a simulation 
experiment can be reported and analyzed is described. The test 
lbed is based on the regression adjusted graphics and estimation 
|methodology developed by Heidelberger and Lewis [Ref. lJ IGG 
lregenerative simulation. From the graphics and associated 
numerics, the experimenter can summarize and see simultaneously 
relative properties, such as bias, normality and standard 







FORM . 
€ 1 § ? 
DD isan 73 1473 = EoITION oF | Nov 4818 oMsoLere UNCLASSIFIED 
/ = e = 
S/N 0102- LF+ O14- 6601 1 secUmity CLASSIFICATION OF THIS PAGE (When Dela Anterec: 





UNCLASD Lr aU 


a I Pp a i 
SECURITY CLASSIFICATION OF THIS PAGE (When Data Entered) 


w20)- ABSTRACT — (CONTINUED) 





deviation, of several estimators of a characteristic of 

a population for up to § sample sizes. The evolution of 
these properties with sample size is also displayed. 

The graphics 1s supported on a line printer to make it and 
the program portable. The technique is illustrated by 
examples concerning the effects of changes in data distri- 
bution on the behavior of the lag one serial correlation 
coefficient, the estimation of the shape parameter of Gamma 
random variables and a comparison of different methods 
(jackknife, bootstrap) for estimating the standard error of 
an estimator. 


S-N 0102- LF- 014-6601 UNCLASSIFIED 


gf SS SSS 
SECURITY CLASSIFICATION OF THIS PAGE(When Data Entered) 





Approved for public release; distribution unlimited. 


SIMTBED a Graphical Test Bed for Analyzing 
and Reporting the Results of a 
Statistical Simulation Experiment 


by 
Hans-Walter Drueg 


Major, German Air Force 
Ing.(grad.), Hochschule der Bundeswehr Muenchen, 1974 


Submitted in partial fulfillment of the 
requirements for the degree of 


MASTER OF SCIENCE IN OPERATIONS RESEARCH 


from the 


NAVAL POSTGRADUATE SCHOOL 
September 1983 





ABSTRACT 


A graphical test bed in which the results of a simula- 
tion experiment can be reported and analyzed is described. 
The test bed is based on the regression adjusted graphics 
and estimation methodology developed by Heidelberger and 
Lewis [Ref. 1] for regenerative simulation. From the 
graphics and associated numerics, the experimenter can sum- 
marize and see simultaneously relative properties, such as 
bias, normality and standard deviation, of several estimators 
mea characteristic of a population for up to 8 sample sizes. 
The evolution of these properties with sample size is also 
displayed. The graphics is supported on a line printer to 
make it and the program portable. The technique is illus- 
trated by examples concerning the effects of changes in 
data distribution on the behavior of the lag one serial 
correlation coefficient, the estimation of the shape 
parameter of Gamma random variables and a comparison of 
different methods (jackknife, bootstrap) for estimating the 


Standard error of an estimator. 
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i eo NOPSZS 


SIMTBED 


THE PROGRAM: 
Portable FORTRAN program using printer plot graphics 
(3 different program versions) 
Program will run on: 
IBM 
VAX 
mer PC 
Scc. 
Ga. 900 lines of FORTRAN Code 
Memory requirements: 
SIMTB1 1 M Bytes 
SIMTB2 1 M Bytes 
SaiMTB 3 0.5 M Bytes 
(may Slightly differ with different type of estimator 
functions and subsample sizes) 
meme OSES : 
- To explore the distribution of a statistical estimator 
- To see how that distribution changes with sample size 
= To compare that distribution with the distribution of 


competing estimators 


ie 





fee USER SUPPLIES: 

A. Optional Parameters: 

NE(1),NE(2),...,NE(8) = Subsample Sizes (maximum is 8) 
The estimator is computed based on NE(i) data points 

N = Total Number of simulated data points per 

replication 

PepouDSam@ple waaze NE(i), there are [N/NE(i)| 
independent values of the estimator 

M = Number of Replications 
When all replications have been run, there are 
M* JN/NE(i)} independent values of the estimator 
at each NE(i) 


D 


Degree of Regression (maximum D = 3) 
L = Number of Subsample Sizes (maximum L = 8} 


GRAPHICS and SCALING options 


Ene Data: 

A total of M*N simulated data values are needed. In 
SIMTB1 and SIMTB2 the same data is used at each 
Subsample size (NE(i)). In SIMTB3 new data is always 


generated. 


fo Estimators: 

Up to 3 FORTRAN functions (i.e. Estimators) are needed. 
They must accept as inputs a data subsample and the size 
of that subsample. They must return one value of the 


estimator for that subsample. 


WA 





Serres PRODUCES: 

- A one page graphical output (box plots) at each 
subsample size 

- Numerical Summaries at each subsample size (mean, 
Std. Dev., Std. Dev. of the mean, skewness, kurtosis) 

= Regressions to quantify changes in the mean and 


variance as subsample size changes 


eS 


4 





ie rine RODUCT EON 


SIMTBED, with the different versions (SIMTB1, SIMTB2 
and SIMTB3) is a graphical display program. The program 
gs based on the program RAGE [Ref. 2]. It is used, with 
Simulated data, on a digital computer. The program can be 
used to examine statistical estimators of different type, 

Or properties of a single estimator under different distri- 
butional assumptions. The distribution of the estimator can 
be explored for given sample sizes and the properties can be 
compared for different sample sizes. The estimation condi- 
Bloons are controlled by the experimenter. It is also possible 
to examine the effects of changes in the underlying distri- 
bution of the data. 

When the program SIMTBi or SIMTB2 is used with simulated 
data, the data is assumed to be independent and identically 
distributed (iid). This iid data can be sectioned into 
M independent blocks of specified sample size N. The sample 
of size N, will be sectioned into smaller subsample of size 
NE(k). The estimates are then calculated from this subsample 
of size NE(k). 

One salient feature of the program versions SIMTBI1 and 
SIMTB2 is that they use the same batch of simulated random 
Variables to explore the properties of all the estimators 


at the various subsample sizes. This is done for economy 


14 





fmemcOuld be important on slow computers; the price paid is 
Maes the regression analysis provided by SIMTBl1 and SIMTB2 
of its graphical output is performed on correlated samples. 

The version SIMTB3 uses one dimensional data and does 
not have this repetition feature. New data is used for each 
calculation of each estimator at each subsample size. The 
data 1s generated when the estimator function is called and 
only the needed batch of the exact subsample size is gener- 
ated. This technique reduces the memory requirements. 

Moreover all data sets are uncorrelated and a much more 
precise correlation can be performed if required. But this 
technique increases the computer time. 

To use the program it 1S necessary only to define the 
optional parameters (see Section IV), supply the simulated 
random variables or a batch of data points, and provide the 
FORTRAN functions for the calculation of the estimators 
which are of interest. The program itself will subdivide 
the input data and feed the data properly into the functions, 
scale the graphic display, produce boxplots and summary 
statistics. A regression will be computed for the mean and 
variance of each estimator based on inverse subsample size. 
The result of the regression is displayed graphically and 
numerically. 

The program is written in ANSI Standard FORTRAN (X3.9- 
1966) and extensively tested on an IBM 3033 computer using 


FORTRAN IV (H Extended) or FORTRAN IV (Gl) compilers. The 


i 





program SIMTBED (all versions) provides all subroutines 

used inside the program. Besides the estimator functions, 
the user has to provide NO additional subroutines via 
packages like IMSL. The program should be totally portable; 
it has been tested under FORTRAN 77 on a VAX 11/780. The 


only limitation is the available memory. 
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III. GENERAL IDEA AND DATA STRUCTURE OF THE PROGRAM 


The main purpose of SIMTBED, with the different versions, 
is to explore the distributional behavior of estimators 


and show their properties in a graphical and numerical 
display. All versions use the same ideas, they only differ 


in the type of data and the way the data is used inside the 





program. 
_ ae 
! 
. 
M 


k 2123, 51 


Figure 1. Sectioning of the Data into M*N Sections 
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To study the behavior of an estimator the experimenter 
usually uses well known simulated data. Thus if one is 
interested in exploring the behavior of estimates of the 
shape parameter in a Gamma population, one generates Gamma 
variates from a random number generator package (e.g., IMSL 
subroutine Chap. G). This batch of simulated data is 
processed by SIMTBED in the following way. 

The data batch is first divided into M independent 
blocks. Each block contains N data points. So the starting 
data batch has to consist of M*N data points (see Figure 1). 

All blocks are divided into subsamples of size n.. The 
actual subsample size n. is an element of the subsample 
Size array NE. This array can store up to 8 different 
values. Then the estimator is calculated for each subsample 
of size n,. The estimator function will be calculated 
([N/n; J) *™ times. This total population of estimates is used 
to evaluate the summary statistics for the estimate and 
construct the corresponding box plot. If NE contains 
another element, the blocks are divided into the new sub- 


sample n. size and all calculations are done again. 


vet 1 
In addition to the summary statistics and the box 
plots (see e.g., Figure 3a), a regression on the averages 
and on the variance is computed, following the methodology 


of regression adjusted estimate (RARE) developed by 


Heidelberger and Lewis [Ref. 1]. 


ie 





The RARE estimate is the regression coefficient Che « 
It is the asymptotic estimate of the expected value of the 
parameter. The unbiased RAGE estimate of the average of the 


parameter is determined by the regression formula: 


7 l l l 
E(9 (n; ) ) = 1K + i n. + m5 > Tr oeecee. oe On an) 


where D is an input parameter. 

The asymptotic RARE estimate is printed as a dashed line 
in the graph. 

The unbiased RAGE estimate of the variance of the esti- 


mator is given by the formula: 


2 _ i 1 m 
S (n;) = By n. t B 5 [aes oe B 5 ae See ia ete 
a (n.) (n.) 


@Mhe regression calculations of the average are done for 
each block separately. The result are M sets of regression 
coefficients. The finally printed regression coefficients 
are the averages of the M replications (see Figure 2). From 
the set of coefficients the variance and the standard 
deviation of the regression coefficients is calculated. 

The regression on the variance is done once, using all 


(M*N) data points. 


Le 








Ware vaRial- | 


OINE(L I 


Regression on 


Figure 2. Structure of the Regression Coefficients 
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IV. ARGUMENTS OF THE PROGRAMS SIMTB1 AND SIMTB2 


All versions of SIMTBED have in general the same argu- 
ment list. SIMTB2 (a version for multivariate random varia- 
bles) has two additional arguments. All arguments will be 
meserioecd in detail and all restrictions or limitations will 
be discussed. In Section V, detailed examples of setups 


Pe oLMTBl and SIMTB2 are given. 


A. X=--DATA ARRAY 

The array X is the input data array. xX is single pre- 
cision REAL*4, and is generally simulated data e.g., in 
Section VII, independent Gamma variates. 

maesiMTBl the dimension of the array X must not exceed 
207000. 

In SIMTB2 the input data is multivariate and conse- 
quently the array X has two dimensions. The size of the 
dimensions is not directly limited by the program. The 
space must be provided by the calling program and is passed 
aS an argument (see IR and IRK). The memory requirements 


increase rapidly as the dimensions of X increase. 


eae N—-—-SAMPLE SIZE 

The sample size N is the number of data points per section 
Peminput data X. N is an INTEGER. Depending on the pre- 
cision of the simulation that is required, the sample size 


Meean vary from 1 to 50,000. 
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In SIMTB2 N is the number of multivariate data points per 
section. M (the number of replications) times N must not 
exceed IR, the row dimension of X. 

If M times N exceeds 50,000 an error message will be 
Peemced and the execution will be terminated. If the product 
M times N exceeds the total number of data points provided 
by the user, NO error message will signal the user error and 


the result of the execution is not predictable. 


fee i—-—- NUMBER OF REPLICATIONS 

The number of replications M of the array X is an 
INTEGER. M determines the number of sections, into which 
the data set X is divided. So M*N is the dimension of X in 
SIMTB1. The parameter M also determines the number of 
regressions on the mean values that will be run to find the 
regression coefficients in the regression on the mean. If 
Mis l only one regression will be done and no variance and 
Standard deviation of the regression coefficients can be 


calculated. 


Be WNE--SUBSAMPLE SIZE ARRAY 

The argument NE is an INTEGER array, containing up to 
8 subsample sizes. These are the subsample sizes at which 
the properties of the estimator are to be investigated. NE 
must always contain 8 elements. If less than 38 subsample 
Sizes are used, the array must be filled up with dummy 


values. The estimator is calculated for each subsample 
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individually. The used values of NE (not the dummy values) 
have to be ordered from the smallest to the largest 
Mme) 6 6 < 6 NE(2) < NE(3) < ...). 

The program can only handle up to a total number of 
12,500 estimates. This limit is caused by storage limita- 
tions and can be expanded for larger computers. The smallest 
subsample size NE(1l) produces the most estimates and if this 
case exceeds 12,500 the execution will be terminated and an 
error message will be printed. 


In SIMTB2 all values of NE have to be less than 5,000. 


E. L--~NUMBER OF SUBSAMPLE SIZES (BOX PLOTS) 

The number of box plots L in the graphical output (the 
number of subsamples at which the estimator is examined), 
is an INTEGER with values from 1 to 8. The value of L 
determines how many elements of the array NE the program 
will use (e.g., L = 2, the program uses only NE(1) and NE(2) 
for the calculations). L determines the number of box plots 
and the number of corresponding summary statistics in the 
output. If L is out of range the program will terminate 
execution and print an error message. 

For L = 1 no regression is possible. No regression 


@eeput will be printed. 


mee D--DEGREE OF REGRESSION 
The degree of regression on the mean D, is an INTEGER 


with values from 1 to 3. The chosen degree refers to the 
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member of coefficients calculated and printed for the 
Begression equations. Experience has shown that D = 3 is 
preferable (higher values cause severe numerical problems in 
the regression computations). 

If D exceeds the value of L-l1, the program reduces D 


to this number, regardless of the value chosen by the user. 


eee nG-~-REDUCED GRAPHICS 
The argument RG is an INTEGER with value: 
MeeeoLt ==> full graphics 
1 = set ==> reduced graphics 
In the reduced graphics case the vertical scale of the 
graph is reduced by eliminating the outliers. Only the 
number of outliers is counted and the number is printed. 
This enables the user to see the body of the boxplot in more 


detail. 


peo L-—-SCALING ESTIMATORS INDIVIDUALLY 

The argument SEI is an INTEGER with value: 
0 = off ==> common scaling 
l = set ==> individual scaling. 

For SEI = 0 all printed graphs (max. 3 per program 
run) are scaled to the same range. This makes the compari- 
son of the different graphs easier. The value SEI = 1 causes 
the program to scale each graph individually. 

The argument SEI works in combination with the arguments 
RG and SVS. The combination of these three arguments make it 


possible to fit the printed graphs to the needs of the user. 
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eo vo--SoLLTING THE VERTICAL SCALE 

The argument SVS is an INTEGER with value: 
0 = off ==> automatic scaling 
1 = set ==> scaling by the user. 

SVS = 0 causes the program to calculate the scale of 
e@emdreaphic printout. The final graphic display is influ- 
enced by the chosen values of RG and SEI. The values of YMIN 
and YMAX are ignored by the program. 

For SVS = 1 the user must provide the scaling. The 
graphics are vertically scaled to a given minimum value 
(YMIN) and maximum value (YMIN and YMAX). The arguments of 
RG and SEI are ignored. Only the numbers of outliers outside 


the display are printed for each box plot. 


eee LN—-—-MINIMUM VALUE OF THE VERTICAL SCALE 

The scaling parameter YMIN is data type REAL*4. It 1s 
the lower limit of the vertical scale. It affects the 
scaling only in combination with SVS = l. If the chosen 
value is so large that the value YMIN lies inside the body 
of the box plot, an error occurs and the program ends with 


an abnormal ending. 


K. YMAX--MAXIMUM VALUE OF THE VERTICAL SCALE 

The scaling parameter YMAX is data type REAL*4. It is 
the upper limit of the vertical scale. The scaling is only 
effected if SVS has value l. If the value of YMAX is too 
small and it lies inside the body of the box plot, the 


program comes to an abnormal ending. 
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YMIN and YMAX in combination with SVS = 1 should only 
Memased for well known graphic ouptut. With this option, 
it is possible to scale the output of different program 
runs to a common scale. In particular if more than 3 esti- 
mators have to be estimated and compared, so that a common 


scale is needed, this option may be used. 


Peo L--NUMBER OF ESTIMATORS 

The parameter NEST is an INTEGER with the value l, 2 
or 3. The value of NEST determines the number of different 
one-page graphic displays the program produces with one call 
from the calling program. Usually the value of NEST is 
equal to the number of different estimators used in the 
program. 

In SIMTB2 the same estimator may be used with different 


(e.g., normal, exponential etc.) distributed data sets. 


Peeeroltl, EST2Z, EST3--ESTIMATOR FUNCTIONS 

These 3 parameters are data type REAL*4 and are used 
EO pass the EXTERNAL function names to the program. An 
external declared function is a function which computes the 
value of an estimator. This function may call other routines, 
but the final value of the estimator must be passed over 
this function name. 

For SIMTB1 each function must have the two arguments X 
and NEK (e.g., FUNCTION VARIANCE (X,NEK)). xX is the data 


array and NEK is the number of data elements in X. 
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SIMTB2 needs for each function four arguments X, NEK, 
mek, IRK (e.g., FUNCTION CORRELATION (X,NEK,IDR,IRK)). Xx 
is a two dimensional data array and NEK is the subsample 
size, for which the function is evaluated. IDR and IRK are 
the dimensions of the array X. | 

If the user wants to use less than 3 estimator functions, 
he must use dummy arguments and choose the correct value of 
NEST. The easiest way to do this is to use a function of 
a previous used estimator again (for details of the pro- 


gramming see Section V). 


ert ol, LTTL2, TTL2--DESCRIPTION OF THE ESTIMATORS 

These arrays are used to pass titles from the calling 
eregram to SIMTB. Each title has to be 120 characters long 
(15 fields, each 8 characters wide). The title is printed 
below the output of the ccrresponding function. 

If the user uses less than 3 functions and titles, he 
has tc use dummy arguments (for details see Section V). 

If the titles are not in the correct format, correspond- 
ing to the FORMAT statement the program will not execute 


properly. 


See LR, LRK--DIMENSIONS OF THE ARRAY X 

These arguments are used ONLY in SIMTB2. IR (row dimen- 
Sion) and IRK (column dimension) are INTEGERS and have to 
match the dimension of X in the calling program. 

If the dimensions don't match, NO error message will be 
mBeeduced and the output is unpredictable. The error may 


not be obvious, and the ouput may look reasonable. 
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V. DUMMY EXAMPLES OF IMPLEMENTING SIMTB AND SIMTB2 
INTO A DRIVER-PROGRAM 

The following examples show how the programs SIMTB1 or 
SIMTB2 can be implemented into a FORTRAN driver-program. 
The only purpose of the examples is to clarify the FORTRAN 
implementation, to avoid programming errors by the user. 
Additional comment lines are added to the program examples. 

In the first example SIMTBl1 is used to compare two 
different estimators (VAR and UNBVAR) of the variance of a 
normal sample. The sample may be generated with a random 
number generator (e.g., LLRANDOMIT). 

In the second example SIMTB2 is used to compare the 
distribution of two estimators. The estimators are the 
covariance (COV) and the correlation coefficient (CORR) of 


a bivariate standard normal sample. 


A. EXAMPLE 1 USING SIMTBL 
MAIN 
Seon AMPLE of SIMTB1 Calling program, it has not to be 
the MAIN program 
REAL*4 X(50000), YMIN, YMAX, VAR, UNBVAR 
Bees bl (lS) t2 (15), £S(i5) 


PYTEGER N, M, NE(8), L, D, RG, SEL, SVS, NEST 


DATA N Y 22007 


28 





DATA M 7 ee 
eewaeNE /10,20,25,50,100,0,0,0/ 
Array NE must have 8 elements, if only 5 are used, add 
dummy variables and set L = 5 
DATA L if 57 
DATA RG / O/ 
DATA SEI / Oy 
DATA SVS / Oy, 
DATA NEST/  _2/ 
BatA Ti /'ESTIMATE','OF THE V', 'ARIANCE','USING ', 
peeve =(1/7N)', '*SUM(X(I',')-XBAR) *', '*2 oe oo ee 
BATA T2 /"ESTIMATE','OF THE V','ARIANCE ','USING ‘', 
eae,  —-N) ) *SUM', *(X(I)—-XB','AR)**2 ‘',7*' '/ 
All 15 fields (each 8 characters) have to be 


initialized 
EXTERNAL VAR, UNBVAR 


Generate M*N independent Normal (0,1) Random numbers 


and store into X 


eel (xX,N,M,NE,L,D,RG,SEL,SVS, YMIN, YMAX,NEST 
PAR, cl UNBVAR,T2,VAR,TL) 


EST3=VAR and TTL3=Tl used as dummy variables 


preOP 


END 
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Q 


FUNCTION VAR (X,NEK) 


Function to calculate the Variance. 


mer calculations inside the function should be done in 


DOUBLE PRECISION 


REAL*4 X(N), VAR 


REAL*8 SUM, XBAR, DVAR 


SUM=0 .0D0 


DO 


10 


DO 


20 


10 I=1,N 

SUM=SUM+DBLE (X(T) ) 

CONTINUE 

XBAR=SUM/FLOAT (N) 

SUM=0 .0D0 

20 I=1,N 
SUM=SUM+ ( (DBLE(X(I) ))-XBAR) **2 
CONTINUE 

DVAR=SUM/ FLOAT (N) 


VAR=SNGL ( DVAR) 


TURN 


END 


FUNCTION UNBVAR (X,NEK) 


C Function to calculate the Variance. 


Geena! calculations inside the function should be done in 
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C DOUBLE PRECISION 
c 
REAL*4 X(N), UNBVAR 


REAL*8 SUM, XBAR, DUNVAR 


SUM=0.0D0 
Bo 10 I=1,N 
SUM=SUM+DBLE (X(T) ) 
fo CONTINUE 
XBAR=SUM/FLOAT (N) 
SUM=0 .0D0 
DO 20 I=1,N 
SUM=SUM+ ( (DBLE(X(I)))-XBAR) **2 
20 CONTINUE 
DUNVAR=SUM/FLOAT (N) 


UNBVAR=SNGL (DUNVAR) 


RETURN 


END 


B. EXAMPLE 2 USING SIMTB2 


MAIN 
EXAMPLE of SIMTB2 Calling program, it has not to be 
the MAIN program 
Berean (ZoO00 2), YMIN, YMAX, COV, CORR 


wee «Sen (15), T2(15), T3(15) 


nl 





PimaGceR N, M, NES), L, D, RG, SEI, SVS, NEST, IR, 


DATA 


DATA 


DATA 


DATA 


IRK 
N J 22007 
M ey 


ire 2000/7 


et yy 


IR and IRK must be equal to the dimensions of X 


DATA NE /10,20,25,50,83,100,125,250/ 

DATA L / 8/ 

DATA RG / 0/ 

DATA SEI / 0/ 

DATA SVS / 1/ 

DATA YMIN/ 0.0/ 

DATA YMAX/ 1.0/ 

DATA NEST/ =.2/ 

eA Tt / ESTIMATE’, 'OF THE C', ‘OVARIANC', 'E er 
fae! '/ 

DATA T2 /'ESTIMATE','OF THE C','ORRELAT1','ON COEFF', 
faelereNG 6", '10%' */ 


All 15 fields (each 8 characters) have to be 


initialized 


PERN AL COV, CORR 


Generate M*N pairs of independent random bivariate 


numbers, each pair being independent N(0,1) random 


ow 





Qo Q ™ 


variables and store into X 


SLL SmiM?’s2 (xX,N,M,NE,L,D,RG,SEI,SVS,YMIN, YMAX ,NEXT 
+ ROO LP COnRR, 2, COV, Tl, IR,1 RK) 


EST3=COV and TTL3=Tl used as dummy variables 


= OP 


END 


MBUNCTION COV (X,NEK, IDR, IRK) 
Function to calculate the Covariance. 
All calculations inside the function should be done in 


E@SUBLE PRECISION 


REAL*4 X(IDR,IRK), COV 


REAL~o SoUM),SUM2 ,SUM3 ,XBARI ,XBAR2 ,EXIX2 ,DCOV 


SUM1=0 .QD0 

SUM2=0 .0D0 

SUM3=0 .QDOQ 

EOe0 2=1,N 
SUM1=SUM1+DBLE(X(I,1)) 
SUM2=SUM2+DBLE (X(1I,2)) 
SUM3=SUM3+DBLE (X(I,1) *X(I,2) ) 

10 CONTINUE 


XBAR1L=SUM1/FLOAT (N) 
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10 


XBAR2=SUM2/FLOAT (N) 
EX1X2=SUM3/FLOAT (N) 
DCOV=EX1LX2- (XBAR1*XBARZ2 ) 


COV=SNGL (DCOV) 


RETURN 


END 


MmenGlION CORR (X,NEK, IDR, IRK) 


mumetion to calculate the Correlation coefficient 


Al 


1 calculations inside the function should be done in 


BEUBLE PRECISION 


— 


hia 4 X(IDR,IRK), CORR 
ReAL*s SUML,SUM2,SUM3,XBAR1,XBAR2Z,EX1X2,VAR1,VARZ, 


COV, DCORR 


SUM1=0 .0D0 

SUM2=0 .0D0 

SUM3=0 .0D0 

DO 10 I=1,N 
SUML=SUM1+DBLE(X(I,1)) 
SUM2=SUM2+DBLE(X(I,2) ) 
SUMeeoUMe rpm (x(t, 1) *xX(1,2) ) 

CONTINUE 


XBAR1=SUM1/FLOAT (N) 
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20 


XBAR2-SUM2/FLOAT (N) 
EX1X2=SUM3/FLOAT (N) 
SUM1=0.0DO 
SUM2=0.0D0 
bo 20 I=1,N 
SUM1=SUM1+DBLE (X(I,1) **2) 
SUM2=SUM2+DBLE(X(I,2) **2) 
CONT INUE 
VAR1= (SUM1/FLOAT (N) ) - (XBAR1* *2) 
VAR2= (SUM2/FLOAT (N) ) - (XBAR2**2) 
COV=EX1X2- (XBAR1*XBAR2) 
DCORR=COV/ ( (VAR1*VAR2) **0.5) 


CORR=SNGL (DCORR) 


RETURN 


END 
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fe olUDY OF THE BEHAVIOR OF SERIAL CORRELATION 
BSoUevethoee One ULFRERENT DISTRIBUTIONS 

PeeCALCULATION OF THE FIRST SERIAL CORRELATION COEFFICIENT 

It is known that for an independent sample from a 
population with finite variance, the distribution of the 
serial correlation coefficient (Anderson and Walker, 1964) 
[Ref. 3] is asymptotically Normal with mean zero and vari- 
ances l/n, where n is the sample size. If the population is 
tiit.da Normal then the bias is exactly -1/n. Since those 
asymptotic properties are frequently used as approximations 
in tests of significance, 1t is important to know how valid 
the approximation would be in small samples from a variety 
of distributions. We will look at that question in the next 
two sections and then go on to consider two alternative 
measures of correlation, Fisher's z-transform and the 2-fold 
jackknifed estimate of the correlation. Their ability to 
reduce bias and/or induce Normality will be examined against 
other changes in the distribution of the estimators, particu- 
larly variance inflation. A simulation study, without 
graphics, of some of these problems was conducted by Cox 


me2o6) {Ref. 4]. 


B. SIMTB1L OUTPUT FOR SERIAL CORRELATION 
Figure 3(a) shows the simulated distribution and sample 


properties of the serial correlation coefficient estimate 


36 





Q LL ei Xp) 
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j=l 
where: 
x } 
> = ee aes 
0 ai 5 
7 e 
Xx = Ke/ (n=l), and 
1 4321 7 
ee n 
ii 2 X./(n-1) 


J=2 


for various sub-sample sizes n= n.- This definition matches 


that used by Anderson and Walker (1964). We consider first 
subsamples of size = 10, and then of size n. = 20, Nn. = S07, 
n, = 40, ne = 50, Ne = 7s nz = 100 and Ng = 150, successively. 


For each subsample size the input sample of N = 5000 simulated 
Normal (0,1) random variables is divided into as many full 
subsamples of size n, as possible, and the serial correlation 
is computed for each of the IN/n; | subsamples of size n;- 
The entire procedure is then replicated M = 10 times, each 
time with a new simulated sample of N = 5000 Normal (0,1) 
variables. 

After all M replications have been run, all the estimates 
of serial correlation for each subsample size are grouped 


together and their simulated distribution is presented via a 
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boxplot and summary statistics. The boxplot follows the 
standards discussed in Mosteller and Tukey (1977) [Ref. 5] 
with the median denoted by a + within the box, the mean by a 
* within the box, the outliers by 0's, and the far outliers 
by *'s beyond the whiskers. The summary statistics include 
the sample mean, sample standard deviation, estimated standard 
deviation of the sample mean (i.e., sample standard deviation/ 
sqrt (M |N/n.J), sample skewness and sample kurtosis of the 
correlation estimates. 

Moeoning at the output, the first (leftmost) boxplot in 


the graph in Figure 3(a) shows the distribution of 


(Simulation 
; Sample Size = 5000 
(# Replications) x (Subsample, = 10x Ea 
Size 


10 x 500 = 5000 


estimates of serial correlation from independent subsamples 
of size n, = 10. Summary statistics for the boxplot can be 
found below the graph in the column labeled "Subsample Size 
10," so that the average serial correlation is -.1074, and 
the estimated standard deviation is .2996. The estimated 
Standard deviation of the serial correlation estimate 1s 
.2996/Y(5000) = .00424. Recall that this refers to correla- 
tion estimates based on subsamples of size 10. 


Since the X-axis of the graph represents subsample size, 


the last (rightmost) boxplot shows the distribution of 
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10 x ao =n) x 33. = sa0 


estimates of serial correlation from independent subsamples 
of size Ng = 150. Although the 330 estimates are independent 
of each other, they are not independent of the 5000 estimates 
that comprise the first boxplot since the same data (divided 
and processed in different ways) was used for both. Summary 
statistics show that the average correlation has dropped to 
-.007372, indicating the fall off in bias, and the standard 
deviation has dropped to .07822, indicating the greater 
accuracy with which the correlation can be estimated when 150 
points, rather than 10, are available. 

In order to quantify the changes that are occurring in 
the mean and variance of the distribution of the estimator 
as subsample size changes, SIMTBl performs two types of 
regressions. The first regression is on the averages and is 
done after each replication, using the average serial corre- 
faeron for that replication, roo as the dependent variable. 
Inverse powers of the Saneeaae sive serve as the indepen- 
dent variables. For Figure 3(a) the degree of the regression 


was chosen to be D = 3 so, for each replication, the equations 


we attempt to fit by least squares are: 


a a a 
= = 1 2 5 here 
Gon a a0 ene ca = Tor ep == We ue er cey Ot 
8 a atc Os 


This form anticipates the general asymptotic expansion 


a 





a O O OL 
me ()) = 6+ —+ ++ + 
n 2 3 
n n 
maeen holds true in the current situation with 90 = 0 and 
(in the Normal case) a = -1l1 (see Cramer (1948) for general 


results of this type) [Ref. 6]. 


Values of Bas Ayr a5 and a, are calculated after each 


replication, averaged across the M replications to get age 


aye as. and as, and then the averages are reported below the 


Summary statistics on the line "Mean of Regression on Averages-- 


Meerticients." We find that ay = —,000272 and ay =o 74 


both close to their theoretical counterparts. 
Because we have 10 replications and therefore 10 inde- 


pendent values of each of a a anda we can also 


Oe mao Bn 


estimate the variances and standard deviations of a a a 


Opa 


and a, across replications. These values are presented on 
the two lines immediately below the coefficients. For 
instance, the estimated s.d. of the estimate ay = -.000272 


or Ao om. 003892). 

The regression line for the mean value of the estimator 
is presented visually in the graph as a dotted curve. The 
estimated asymptote (i.e., ay) is printed with a dashed line 
wherever it does not coincide with the regression line. Bias, 
therefore, can be viewed as the difference between those two 
lines. 


The second regression referred to above is done after all 


replications have been run and the variances of the estimators 
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at each subsample size have been calculated. (Note that the 
standard deviations, not the variances, are presented in the 
summary statistics.) It should be recalled from previous 
discussion that these variances, as well as all measures in 
the summary statistics, are based on the grouping together 
fee Serial correlations from all replications, at each 
subsample size. This 1S in contrast to the procedure for the 
regression on the means, where average correlations are com- 
Meee cOr Gach subsample size for each replication. In the 


case of the variances, we have 8 equations: 


ee Oy 


-- 
lI 








which we fit by least squares in order to estimate the coeffi- 
cient Bye By: Bos and 8, in the presumed asymptotic expansion 
Bg By B 5 B 3 


Var(O(n)) = ot 372 + a 572 ee 








This expansion holds for the variance of the estimated 
serial correlation coefficient for independent data. Usually 
it will be By in which we are most interested since By is 
used in computing asymptotic relative efficiencies of esti- 
mators. For independent data with finite variance, we Know 
that Bo = 1. The computed values of Doe boy bo, and bo, 
are presented on the line labeled ‘Regression on Variance-- 
Coefficients'. Notice that b, = .7438 is close to the 


0 


theoretical value of l. 
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The final two numbers on Figure 3(a), YMIN and YMAXx, 
Simply show the scale of the vertical axis. Because the SIMTBL 
program option to put Figures 3(a), 3(b) and 3(c) on the 
same scale was in effect, it may be that no boxplot in a 
@vyen Figure (e.g., Figure 3(b)) requires the full range of 
Y-values. 

In order to produce Figure 3(b), the Normal (0,1) data 
that went into Figure 1(a) was squared to create longer 
tailed (i) random variables. The output is entirely analo- 
mememco that for Figure 3(a). Similarly, for Figure l(c), 
the Normal (0,1) data was exponentiated in order to create 
Lognormal (0,1) data and to produce analogous graphical 
output. The indication is that the distribution of the sample 
serial correlation is robust with respect to the population 
Seoesabution. 

The features of the SIMTB1 output will become clearer when 
they are associated with the various properties of the corre- 
lation estimator. First, however, a few technical comments 


concerning the regressions are necessary. 


C. SOME COMMENTS ON THE REGRESSIONS 

Two types of problems, numerical and statistical, can 
occur when attempting to fit the two sets of regression 
equations presented in Section VI.B. 

First, there is the question of numerical stability when 


the independent variables, {l, — nee n>} 


a ae 
3/2 2 a 


or 


ier ; } decrease geometrically. If we 


— iL ree: 


42 





attempt to form xix, where X is the respective design matrix 


i. 
and X is the transpose of X, we get values that range from 
8 


8 (assuming 8 subsample sizes) to } ae for the regression 
8 8 1=1 
on the means, and } ni comme i for the regression on 
1=1 1=1 


the variances. Experience has shown that attempts to solve 
systems with such extremes in the xt x matrix produce errone- 
ous results. Instead, SIMTB1 scales the design matrices by 
multiplying each column of X by Max(n,) raised to the proper 
power so that no entry becomes nee small. The standard 
Choleski decomposition is then used to fit the equations, 
and the coefficients are properly rescaled before they are 
reported. This procedure produces numerically reliable 
results. 

The second problem concerns the breakdown of statistical 
assumptions in our regression models. It has already been 
pointed out in Section VI.B that the two sets of dependent 
variables: 


(1) the O(n,) when considering the regression on the means; 


5 =~ 2 
Lean a CO tm.) = O(n )) 
(2) the s*(n.) = MS i a = i 
aa iF il 


when considering the regression on the variances, 
are not independent over i since all are based on the same 
Simulated data. The extent of the dependence is demonstrated 
by the correlation matrix in Table l. Entries in that table 


; Z 
show the estimated correlation between s*(n,) and s aoe EO 
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TABLE 1 


ENTRIES IN THE TABLE ARE THE ESTIMATED CORRELATIONS 
BETWEEN THE ESTIMATED VARIANCES OF THE a AT 
DIFFERENT SUBSAMPLE SIZES: a 


CORR(s* (x ee (x )) tes” Tae l,..cgee 4) = 1,8 


a a Z 5 d S 6 7 8 

J 

1 1.00 ~ 49 Ome 226 eS ae ey, 14 al 
2 -49 1.00 . 40 yO) eel 38 75 Cmee Ue 
3 . 46 -40 1.00 125 23 44 al a2o 
4 ae © 2 523. eaesUU ~42 . 86 a7 5 
5 . eS MIL 2S <icme 1.00 7 ~43 209 
6 erly . 38 44 . 86 ee ee 0 nao so35 
y 14 136 ak oe 43 “40° - 1.00 cane 
8 Oe es 29 35 ooo 7o5 pee) 


Recall that ry is the estimated serial correlation for a 

Simulated Normal (0,1) subsample of size n. Also, the 

estimated correlations shown above were computed using 10 

values (replications) of See ) and s*(r_ ) for each i and 
1 J 

al. 
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all i and j, where the estimation was done by repeating the 
SIMTB1 experiment with 10 different batches of 50,000 simu- 
lated random variables. Since only 10 values went into each 
correlation calculation, the table is only accurate to 
within approximately + 2/Y10 = .632. We see some indication 
of positive correlation, especially when i and j are close, 
but the lack of independence is not severe enough to hurt 
the regression results for either the estimated means or 
variances Significantly. 

A second assumption, implicit in any regression, is that 
the dependent variables have equal variances. This condition 


holds true for the means, which can be shown to satisfy 
= M 
Var (O(n,)) = ie 


independently of 1. The estimated variances, however, are not 
equivalent and, if we assume the caer to be approximately 
Normally distributed so that 
Miya] a | 5 
j=l 
’ Z ; 
1S approximately proportional toa XM /ng}- 1 random variable, 


we can compute 


Ie 


ee ee 5 
MNn.-n.: 
ot! 
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To correct this problem, SIMTBl1 scales the s*(n,) by yn, 


So that 


tl 2 


2 
Var (/n, S (n,)) 


Since MN >> n.- The design matrix is scaled accordingly and 


b b and b. discussed in Section VI.B. are 


the values b 5! 3 


Od: 
Beported. 

Table 2 shows the effects of the rescaling by presenting 
first the estimated variances of the s*(n.), where the esti- 
mation is done by repeating SIMTB1 for 10 batches of 50,000 
Simulated data points. These estimated variances decrease 
as n; increases, Closely paralleling the second line of Table 
2 which has the approximate theoretical values (l.e., 
2/(MNn, -n‘)). The final line of Table 2 shows the estimated 
variances of the rescaled s“(n.), i.e., the vn. s*(n.), 
which, as expected and hoped, show a more constant variance 
moth 1. 

Although future versions of SIMTB1 will include more 
sophisticated regression routines and the ability to gener- 


ate independent samples at each subsample size, the SIMTB1 


is quick, usable, and accurate for most situations. 


D. INTERPRETING THE SERIAL CORRELATION RESULTS 

Returning to Figure 3(a) which shows the simulated dis- 
mepoution of the serial correlation coefficient from indepen- 
dent, Normal (0,1) data, the following comments summarize the 


most striking features: 


AL 





TABLE 2 


A COMPARISON OF THE ESTIMATED VARIANCE OF s*(r ) 


WITH THE APPROXIMATE THEORETICAL VARIANCE OF “"j 
s“(r_) AND WITH THE APPROXIMATELY EQUIVARIANT SCALED 
i VERSIONS, the Seam a 


a 
All entries have been multiplied by LO? 


n. = LO 20 30 40 30 7 J2hO0.. 150 


Var(s*(r_ )) iio OmercO4 .07/9 995047 .031 .049 2.022 


Approx. Theoretical 


Var (s*(r_ )) BAO C200 eal oe LO Oi 080 9053, ..040. 027 
of 


2 se P ieee omni ee eto? oo 2. 33k 88. 38.35 


Var(n.° 
‘i 
i 


The estimated variances of ce ) and vn; s*(r_ ) were 
a a 
calculated using 10 independent replications of s*(r_ ive 
1 
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(a) The boxplots appear very symmetric at all subsample 
sizes with nearly equal numbers of outliers at either tail 
and with mean and median coincidental. This observation is 
confirmed by the estimates of skewness in the summary sta- 
tistics. Kurtosis is mildly negative at small subsample 
sizes but, overall, asymptotic Normality seems to take hold 
ieeemer quickly. 

(b) The average serial correlation is negative for small 
subsam-les. This is demonstrated by the dotted regression 
curve which starts at approximately -.10 and levels off near 
0 for subsamples greater than about 85. The dashed asymptote 
of -.000272 is very close to the theoretical value of 0, and 
the mean values in the summary table closely reflect the bias 
of -l/n. 

(c) The standard deviations of the simulated distributions 
are very close to the asymptotic values of eo although 
the lead coefficient in the regression on the variances, 

Do = .743756, is not as close to the theoretical value of l 
as we would hope. When SIMTB1 is repeated 10 times with 10 
different batches of simualted data, we find an average value 
aor by to be 1.0604, with a standard deviation for Bo of 
307. The estimation procedure for bos therefore, remains 
valid, but the estimate itself is highly variable. 

The agreement between the simulated and the theoretical, 


asymptotic values of the bias and variances was discovered 


previously by Cox (1966). SIMTB1 has now allowed us to 
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automatically look at a broader range of subsample sizes and 
to see, through boxplots and estima es of skewness and 
Kurtosis, a fuller picture of any changes in the distribution 
of the estimator. We can be satisfied that estimates of 
serial correlation do behave approximately as Normal (-l/n, 
1/n) rangome variables when the underlying data is Normal 
mor, i). 

If the lead terms in the expansions of the mean and 
variance of the estimated correlation coefficient (i.e., Ags 
ays and bo) had been unknown, we would also have a fairly 
good idea now of what they were. 

When the underlying data is ae Pigiecm so) COmf iris 
Cox's observation that the bias is relatively uneffected but, 
for small subsamples, the standard deviation is smaller than 


ee mime Figure 3(a), Enere asa pronounced 


the expected n_ 
skewness in the boxplots in Figure 3(b) with many more out- 
liers at the positive end, and with the mean higher than the 
Meaian at the first four subsample sizes. The problem of 
Suppressed variance seems cured at no = 100 and Ng = is 0, 
but the skewness remains and could cause problems in tests 
of significance. 

Figure 3(c), which is based on an underlying batch of 
Simulated Lognormal (0,1) data, shows a slight exaggeration 
of the effects in Figure 3(b). The standard deviation is 
more suppressed and does not attain the theoretical level 
by n, = 150. The positive skewness is more pronounced and 


8 


kurtosis does not approach the theoretical value of 0. 
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Overall, the effects of long-tailed data on the distribu- 
mom OL the serial correlation coefficient can be summarized 
as follows: 

(i) Bias 1s not significantly effected and remains at 
approximately -l/n. 

1. ) The variance of the distribution of the serial 
correlation coefficient is reduced by longer- 
talled data. 


(i111) Positive skewness is created in the distribution. 


(av) Kurtosis may become positive at larger subsample 
sizes. 
(v) For long-tailed data (1.e., Lognormal), a subsample 


size of 150 is not large enough to insure 


asymptotic Normality. 


rs. SIMTB1L OUTPUT FOR THE Z-TRANSFORM OF THE CORRELATION 
Fisher's z-transform of the estimated correlated coeffi- 
cient is defined by: 


1 l+e 
C= ae NG 
Z -=- yy. * 
n 


where r, is the estimated serial correlation presented in 
Section VI.B. The transformation is intended to make the 
distribution of the Z more Normal than that of the ro: 
When the same SIMTB1 experiment described in Section VI.B. 
1s run using Z, as the estimator instead of ros we get the 


results shown in Figures 4(a), 4(b) and 4(c). It should 


50 





be noted that the scale of the boxplots here has been forced 
to be approximately comparable ta the scale for the boxplots 
in Section VI.B. This is done by suppressing outliers that 
are more than 1.5 interquartile distances beyond the quar- 
tiles of the boxplot. If we had allowed the data to scale 
the boxplots, we would have seen a much wider range on the 
vertical axis because the Z 4 are not restricted to the limits 
of ~l1 to +l and because there is one far outlier at -3.8. 
In this type of "reduced graphics," we still see the number 
of outliers that fall beyond the allowable range through the 
numbers at the ends of the boxplots, but we do not see their 
actual locations. 

Figure 4(a) shows the distribution of the z-transformed 
correlation coefficients when the underlying data is simu- 
lated, Normal (0,1). At each subsample size, the mean and 


I and 


Seandard deviation are close to the theoretical n- 
/- respectively. The skewness and kurtosis at subsample 
size n, = 10 are far from the theoretical Normal distribution 
values of 0 and 0, reflecting partly tne one far outlier at 
-~3.8 and partly the negative skew in the remainder of the 


Z 's. For other subsample sizes, there is no strong evi- 
i 


dence to contradict the assumption of approximate Normality. 

The relationship between Figure 4(b) and 4(a) is similar 
to that between 3(b) and 3(a). Figure 4(b), which is based 
on simulated x4 data, shows (a) bias that is the same as for 


the transformed correlations based on Normal data, (b) slightly 
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Suppressed variances, particularly at small subsample sizes 
and (c) positive skewness which persists at large subsample 
Sizes. In addition, there are signs of positive kurtosis 
at small subsample sizes. 

Figure 4(c) is based on Lognormal (0,1) data and shows 
high values of skewness and kurtosis at almost all subsample 
sizes. Approximate Normality seems an unwarranted assumption. 
In fact, the kurtosis is converging very slowly to its 
Beympcotic value of O. 

In general, using the z-transform does not help with 
Normality assumptions, especially when dealing with long- 


tailed distributions. 


PS. SIMTB1 OUTPUT FOR THE 2-FOLD JACKKNIFE OF THE CORRELATION 

The final Figures, 5(a), 5(b) and 5(c), deal with the 2- 
fold jackknife estimate of correlation. Again, the figures 
are reduced graphics with scaling comparable to that of the 
boxplots of Sections VI.D and VI.E. To define the estimator, 
we start with a given subsample of size n, compute the serial 
correlation for the first [n/2] Polntswanea call it r,(n/2), 
compute the serial correlation for the second Ba points and 
earl it r,(n/2) and compute the serial correlation for the 
entire subsample of n points and call it XQ (n). Each compu- 
tation follows the formula in Section VI.B. The three esti- 
mators are then combined to form two pseudo-values, 


r.*(n) = 2x, (n) r, (n/2) 


oe 





and 
(i) 2x, (n) - rj(n/2) , 


and the final jackknife estimator for that subsample is de- 


fined as 


Ene Cx * (1m) 


Although a jackknife estimator may have many favorable 
properties, we are concerned here primarily with its ability 
to remove bias, hopefully without inflating the variances 

of the estimator and/or inducing nonnormality. 

Figure 3(a), based on simulated Normal (0,1) data, shows 
nearly complete removal of bias, even at small subsample 
sizes. _The cost of the bias reduction is reflected in an in- 
crease of nearly 50% in the standard deviation of the corre- 
lation estimate for subsample size 10, and lesser relative 
increases at larger subsample sizes. There is also an indica- 
tion of a positive skew for small subsample sizes, and the 
problem that the jackknife estimator need not fall into the 
-l1 to +l range which is desirable for a correlation coeffi- 
cient estimate. 

When using simulated xs data as in Figure 5(b), or simu- 
lated Lognormal (0,1) data as in Figure 5(c), there is again 


nO problem with bias. Variance inflation, though it exists 


DS 





at small subsample sizes, is not as large as when Normal 
mepee) Gata is used. The distributions of the jackknifed 
correlations show very pronounced positive skews, however, as 
well as positive kurtosis. These two problems are worse for 
the longer-tailed Lognormal data. 

Overall, the jackknife estimator is very successful at 
removing bias but the costs include variance inflation, which 
can be severe at small subsample sizes, plus increased posi- 
tive skewness and kurtosis when the estimates are based on 


data from longer-tailed distributions. 


See COMPARISON OF THE THREE ESTIMATES OF CORRELATION 

For Normal (0,1) data, the distribution of the usual 
correlation coefficient displayed in Figure 3(a) behaves 
very much as theoretical asymptotic calculations would pre- 
dict, even at small subsample sizes. This makes it possible 
to correct for bias in the estimator and to perform tests 
of significance. Use of Fisher's z-transform, as illustrated 
in Figure 4(a) does not seem necessary since it does not 
Significantly improve the approximate Normality of the esti- 
Mator. The jackknife estimator in Figure 5(a) may be valuable 
if a direct, unbiased estimator is needed but the inflated 
variance of the jackknife estimator may limit the usefulness 
of the estimate as well as make any tests of significance 
too conservative. 

When the underlying data comes from a longer-tailed dis- 


tribution, the usual correlation coefficient in Figures 3(b) 
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and 3(c) retains a predictable bias term, although the vari- 
ance of its distribution is slightly depressed and the skew- 
ness and kurtosis becomes positive, even for subsamples as 
large as 150. This means that it is still possible to 
estimate the correlation accurately, but tests of signifi- 
cance fall on shakey assumptions of Normality. The z-transform 
in Figures 4(b) and 4(c) does little to firm up those assump- 
tions and, in some cases, makes the situation worse. As in 
the case of Normal data, the 2-fold jackknifed correlation in 
Figures 5(b) and 5(c) is bias-free but follows a fairly non- 
Normal distribution which would invalidate significance 
testing. 

All of the preceding observations and conclusions flow 
immediately from the nine figures presented so far. Further 
studies could easily be done through SIMTB1, looking at larger 
subsample sizes, correlated data, and alternative marginal 
distributions. For demonstration purposes, though, it is 


better to proceed to our second application. 
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vei. colby (Or TPROBEEMS OF ESTIMATING SHAPE 
PARAMETERS FOR HIGHLY SKEWED DISTRIBUTIONS 

A. ESTIMATING THE SHAPE PARAMETER FOR A GAMMA DISTRIBUTION 

As a second application of SIMTBl1, we will consider a 
meeomem which has received much less statistical attention; 
asymptotic results are summarized in Cox and Lewis (1966, 
Ch. 2) [Ref. 7] and Johnson and Kotz (1970, Ch. 17) [Ref. 
8]. We want to estimate the shape parameter, K, for a Gamma 


distribution, where the Gamma density is given by 


Notice that the mean of this distribution is u, not K/u as 
in some differently parameterized versions of the Gamma den- 
sity. For the data that will be simulated for use in SIMTBI 
mw! use K = 5 and u = 1andK= 0.25 and p= tl. The 
meeer the mean of our estimate is to 5 or 0.25, the better 
(in terms of bias) is our estimation procedure. Other factors 
such as the variance and Normality of the estimator will of 
Course also have influence in the determination of a preferred 
estimator. 

Section VII.B will compare the commonly used maximum 


likelihood estimator to the competing method of moments 
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estimator. Both procedures result in asymptotically Normal 
estimators (Cramer, 1948) but the m.l.e. is usually pre- 
ferred because of its favorable asymptotic relative efficiency 
(Cox and Lewis, 1966) [Ref. 7]. Through SIMTBl1, though, we 
will see that for small subsamples the estimated variances 

of the two estimators of K are not as far apart as asymp- 
totic results lead us to believe. In addition, the bias that 
appears in both estimators is smaller for the moment 
estimator. 

In Section VII.C. we will use a four-fold jackknife of 
both the m.l.e. and moments estimators to successfully remove 
the bias. What is remarkable is that, unlike the jackknifing 
of the serial correlation, there is little or no cost in 
terms of variance inflation and nonnormality for the jack- 
knifed moment estimator. When K = .25, we will see in Section 
VII.D. that the jackknifed m.l.e. dominates the other three 
estimators at all subsample sizes when considering the mean, 


variance, and Normality of the estimator. 


B. MAXIMUM LIKELIHOOD AND MOMENT ESTIMATORS OF K 
Figure 6(a) is very similar in format to the figures 
that have already been presented for the correlation example 
except that: 
(1) The estimator whose distribution is being displayed 
is the maximum likelihood estimator of K, the shape parameter 
of a Gamma(5) population. We denote the estimator, computed 


meom a Simulated subsample of size n, by K(n) and define it 
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Boeme the solution of the equation: 


x A n n 
milltog K(n) -¥(K(n))] = on log ) X./n - ) NO) Guse 
alt fal - 


where the X 5 are the simulated Gamma(5) random variables 
and ¥(-) is the digamma function (Cox and Lewis, 1966). 

(2) The eight subsample sizes which we will be looking at 
are n, = o8u n. = Nn, a LOOT ne = eS ne = liG Gy 


na eae o0 and Ng = 500. We will not see as much detail at 


small subsample sizes but we will see some of the asymptotic 


eee ees hh 


me) 200) effects coming in. 
(3) At each subsample size we will work with M* = 20 
independent replications of N* = 2500 simulated Gamma(5) 


Pemeom variables, instead of the M = 10 replications of 

N = 5000 variables used previously. The total number of 
independent simulated random variables across replications 
remain constant at the program maximum of 50,000. Hence, 
the boxplot at subsample size 50 in Figure 6(a) represents 
the distribution of M*|N*/50] = 1000 estimates of K(50) just 
as the boxplot at subsample size 50 in Figure 3(a) represents 
Peer adistribution of M [N/So| = 1000 estimates of r(n). As 
mmoeas the product, MxN, remains constant, the only effect 
that changing the number of replications has, up to rounding 
ia IN/ng| is to change the results in the regression on the 
averages. By uSing M* = 20 and N* = 2500, SIMTB1 reports 


regression coefficients averaged over 20 replications, but, 
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within each replication, the dependent variables are averages 
Seer just 2500/n,| values of the estimator. 

(4) The boxplots are presented using the reduced graphics 
option. In this option any extreme outliers (i.e., those 
beyond 1.5 interquartile distances) are included as a count 
at the tail of each boxplot. This option was chosen in order 
to give more graphical weight to the body of the distributions 
and the fall-off in the bias. Limited printer resolution 
makes it impossible to show details in the body and the tails 
Gieche distributions if there are many straggling outliers. 

In the case of very extreme outliers, no detail would be seen 
in the body of the boxplot without the reduced graphics option. 

Figure 6(b) looks at the distribution of the moment 
estimator of K, the shape parameter of a Gamma (K) population: 

K(n) = (nel) #7 5 (x, -¥)* , 
i=l 
n 

mere <X = e X./n, n is the subsampie size, and the X, 
are the ii tea Gamma(5) random variables. The SIMTBl 
options and parameters mentioned in (2), (3) and (4) preceding 
are also in effect here. 

The two Figures, 6(a) and 6(b), show a very pronounced 
bias in both estimation procedures, although the moment 
estimator is slightly closer to the unbiased value of 5. 


As expected, the standard deviation of the m.l.e. is lower 


than that of the moment estimator although the relative 
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difference at small subsample sizes, for instance 1.448 versus 
1.482 at n, = 33, may not outweight the increase in bias with 
the m.l.e. At larger subsample sizes, the relative differ- 
ence is close to the theoretical asymptotic relative efficiency 
feo (1.ec., .91 at no = Zou 

Both estimators also show distributions with positive 
skewness and kurtosis that decrease to the asymptotic 0 
levels as subsample size increases. The asymptotics appear 
to take hold more quickly for the moment estimator than for 
the m.l.e. 

In summary, SIMTBl1 shows that the m.l.e. is indeed better 
than the moment estimator in terms of variance, but not as 
good for small sample sizes as asymptotic results would lead 


us to believe. In the other areas of bias and asymptotic 


Normality, the moment estimator would have to be preferred. 


feo -FOLD JACKKNIFED ESTIMATORS OF K 

Figures 6(c) and 6(d) show the distributions of the 4-fold 
jackknife m.l.e. of K and 4-fold jackknifed moment estimator 
of K, respectively. A 4-fold jackknife estimator is similar 
to the 2-fold jackknife estimator described in Section VI.F. 
except that there are 4 pseudo-values that come out of dividing 
each subsample into fourths. More details can be found in 
Mosteller and Tukey (1977) [Ref. 6]. 

The purpose of the jackknife is to remove the conspicuous 
bias observed in Figures 6(a) and 6(b). This goal is seen to 


‘be accomplished in Figure 6(c) and 6(d) and we can also note 
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Smaller values of skewness and kurtosis, indicating a quicker 
approach to asymptotic Normality. The skewness and kurtosis 
of the jackknifed moment estimator are the lowest, at small 
subsample sizes, among all estimators. The variance of the 
jackknifed moment estimator is also only slightly inflated, 
as 1s the variance of the jackknifed m.l.e. 

All told, the jackknifed moment estimator, because of 
its lack of bias, small variance, and low skewness and kurtosis, 
would be the method of choice if estimation of K or signifi- 


cance testing was the goal. 


Pees ULTS FOR K = 0.25 

ioerigures 7(a), 7(b), 7(c) and 7(d) we show similar 
results to those discussed above for the case K = 5.0, but 
Memeo K = 0.25. The fact (Cox and Lewis, 1966, Ch. 3) [Ref. 
7] that the m.l.e. estimate is much more efficient than the 
moment estimate is graphically illustrated. What is new is 
the effect of jackknifing: bias 1s reduced without the 
sacrifice of variance inflation or nonnormality. 

Further comparisons and interpretations are similar to 


those done for the case K = 5.0, and are left to the reader. 


ee CONCLUSIONS 

Simply by providing SIMTB1 with the desired estimators, 
we have been able (a) to explore in depth the effects of 
changes in data distribution and of different estimation 


procedures on the calculation of the serial correlation 
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coefficient, and (b) to compare four different ways to esti- 
mate the shape parameter in a highly skewed Gamma population. 
The graphics and numerical output combine to let us see 
and quantify distributional changes that occur as subsample 
size grows. We can see bias fall away, variance shrink, and 
skewness disappear as the estimator approaches asymptotic 
Normality. Terms in the asymptotic expansion of the mean 
and variance of the estimator are automatically calculated 


and can be used to compare different estimators. 
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VIftft. COMPARISON OF DIFFERENT METHODS FOR ESTIMATING 
THE VARIABILITY OF THE STANDARD DEVIATION OF 
CORRELATION ESTIMATES 
A. INTRODUCTION 

Bradley Efron and Gail Gong (1982) review in their 
a@eeacie, "A Leisurely Look at the Bootstrap, the Jackknife, 
and Cross-Validation" [Ref. 9] different methods for esti- 
mating statistical error of parameter estimates. They dis- 
cuss in particular the problem of estimating the error of a 
statistical estimator with the example of the estimates of 
the standard error for the correlation coefficient of a 
Divariate normal distribution. They compared the standard 
deviation estimates of different methods using 200 simulations 
at one fixed sample size of 14. 

SIMTB2 was used to explore the distributions of the 
estimates they used in their article. The bootstrap, the 
jackknife, the infinitesimal jackknife (delta method) and 
the normal theory are the methods for which the distributions 
were explored. The estimation methods will not be explained 
in detail (see Efron and Gong in The American Statistician, 
fees 1983, Vol. 37, No. 1), only the setup of SIMTB2 and 


the output will be discussed. 


Eee sslUP OF SIMTB2 AND THE ESTIMATOR FUNCTIONS 
The bivariate normal distributed input data with known 


correlation 0.5 was generated with an IMSL random number 
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generator (GGNSM). 14,000 bivariate data points were gener- 
ated and then sectioned into 10 blocks, each having 1,400 
points (M = 10, N = 1400). The subsample sizes (NE(k)) used 
ere 10, 14, 20, 28, 35, 40, 70 and 100. Only subsample size 
14 was used by Efron and Gong. 

The program had to be run with 5 different standard 
deviation estimator functions. SIMTB2, as the other ver- 
Sions, can only handle up to 3 estimator functions in one 
program run. To make the given outputs comparable the fixed 
scale option (SVS = 1, YMIN = 0.0 and YMAX = 0.7) was chosen. 

The bootstrap estimate of standard deviation was done with 
2 different numbers of bootstrap replications (B = 128, 

B = 512). So for the subsample size of 14, the bootstrap 
procedure with B bootstrap replications in itself was done 
meuo714 = 100; 100*10 = 1000) 1,000 times. 

The jackknife function followed the standard jackknife 


procedure and used the jackknife formula for the standard 


meVlation: 
n 7 2 
n-l ~ Z 
S. = [—— ) (X,., - X,,,)7] 
jack n 324 a) (.2.) 
with 
24) = O(K) Xone ee Xs ye Kag pres Xp) iO 1e X05) 
and 
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The delta method (infinitesimal jackknife) function 


calculated the estimate with the formula: 





a2 a ws nN “~~ a 
, ee ee 04, “Mame 22 | 32 3_ |”? 
Delta on a2 oa 2 “A “A A ~~ me i 
eae?  Maorgo an aialt eee tet 
with 
A ho h 
4 =\g 
' = = ) (¥, -¥)7(2, -2) 
a nN j=1 
and 


In doing the actual calculations for the delta method, 
for small subsample sizes (10,14,...,40) the value of the 
Variance becomes negative. A negative variance can not be 
interpreted with much meaning. The number of negative values 
goes down (from 305 for NE(1) = 19 to 18 for NE(6) = 40) 
with increasing subsample size. Most of the negative values 
are small. To solve the programming problem (Square root 
of a negative number) the function sets this negative value 
ee 0.0D0. 

For the normal theory estimate, instead of carrying out 


a bootstrap, an approximation formula was used. It is the 
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same formula the article uses for the comparison. 


PE egecs) 
norm i= 3) 3) 


See Efron and Gong for details. 


SeeeeeNLERPRETING THE SIMTB2 OUTPUT 

The output of the program runs is provided as Figures 6, 
fee 2 and 10. A comparison of the numerical output (sub- 
sample size 14) with the results is done in Table 3. 

The bootstrap procedure (Figures 8 and 9) was done with 
2 different numbers of bootstrap replications (B = 128 
and B = 512). Both distributions for the standard deviation 
(S.D.) look very similar. Both are positively skewed with 
some outliers at the right tail. In both cases the outliers 
are in the same range. As the boxplots and the summary 
statistics show, the increase of the number of bootstrap 
replications (B) does not result ina large improvement in 
the performance of the estimation function. 

The jackknife estimate (Figure 10) has a positively skewed 
distribution with outliers. The distribution of the jack- 
knife estimate looks very similar to the distribution of the 
bootstrap estimates. For small subsample sizes the bootstrap 
distribution has more outliers. Overall the performance of 
the jackknife procedure is as good as the bootstrap, but the 


jackknife needs less computer time. 
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The distribution of the estimates produced by the delta 
method (Figure 11) is negatively skewed and has nearly no 
outliers. But in calculating the estimates the problem of 
negative values for the variance came up. For some sub- 
samples, the final estimate (the standard deviation) could 


not be calculated, since the corresponding value of the 


variance was negative. In these cases, the standard devia- 
Puemewas Set to 0.0. This procedure influences the distribu- 
tion and the summary statistics. The influence is more 


important for small subsample sizes than for larger ones. 
So the graphical and numerical output should be seen with 
Peo tact in mind. 

The normal theory function (Figure 12) produces esti- 
mates with a negatively skewed distribution but only a few 
outliers and the distribution is tailed to the left. The 
@oeeeor the distribution is in the opposite direction of 
all other distributions. For the estimate of the standard 
deviation for the correlation coefficient the result of the 
normal theory is close to the true result. This may not be 
valid for other estimators. 

In addition to the comparisons Efron and Gong did, with 
SIMTB2 it is easy to investigate how the sample size will 
influence the estimate of the standard deviation of the corre- 
lation coefficient. In Table 4 the methods are compared for 
a subsample size of 10 and 100. With increasing subsample 


size the quality of the estimate should improve, but the 
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improvement may be different for the different methods of 
estimation. By making the subsample size 10 times larger, 
with the SIMTB2-side effect of reducing the total number of 
calculated estimates, the bootstrap improves less than the 


@eeennire, delta method and normal theory. 


gE8 





TX. FINAL CONCLUSIONS 


SIMTBED with the different versions can be used on 
digital computers of different size (mainframe to micro) 
and type. The limitations in using the program are given 
more by hardware constraints like memory size and computer 
time than by the program itself. 

The FORTRAN program is completely portable, changes in 
the code, may only be necessary to adapt the program to 
special restrictions given by a special type of hardware. 
This may occur in using micro computers more often than 
with mainframe computers. Up to now all versions of code 
are written for the more normal standard computer environment 
and do not need special equipment (color printer, etc.). ° 
Additionally hardware dependent features like color output 
can improve the graphics of the program. 

SIMTBED makes it easy to evaluate the result of statis- 
tical experiments. The combination of graphics and numerical 
summaries for different sample sizes make it easy to judge 
the distributional behavior of a statistical estimator. The 
result can be seen without additional computations in the 
graphic outputs. Comparing only the boxplots it is possible 
to judge the influence of subsample size on the variability 


of an estimator. 


wk 





Besides for research the program can be used in showing 
students the distributional behavior of different estimators 
in a pictorial way. It is easy to compare the different 
behavior of similar estimators (e.g., biased vs. unbiased 
estimator of the variance) for different sample sizes. 

The easy use and the fast visual impression of the 
G@estributional behavior of an estimator, given by the 
Graphic output is one of the advantages in using SIMTBED. 
Besides this fast first visual impression all necessary and 
needed numerics are given for further and deeper 


investigations. 


ee 
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